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ABSTRACT 

Most cool core clusters of galaxies possess active galactic nuclei (AGN) in their centers. 
These AGN inflate buoyant bubbles containing non-thermal radio emitting particles. If 
such bubbles efficiently confine cosmic rays (CR) then this could explain "radio ghosts" 
seen far from cluster centers. We simulate the diffusion of cosmic rays from buoyant 
bubbles inflated by AGN. Our simulations include the effects of the anisotropic parti- 
cle diffusion introduced by magnetic fields. Our models are consistent with the X-ray 
morphology of AGN bubbles, with disruption being suppressed by the magnetic drap- 
ing effect. We conclude that for such magnetic field topologies, a substantial fraction 
of cosmic rays can be confined inside the bubbles on buoyant rise timescales even when 
the parallel diffusivity coefficient is very large. For isotropic diffusion at a comparable 
level, cosmic rays would leak out of the bubbles too rapidly to be consistent with 
radio observations. Thus, the long confinement times associated with the magnetic 
suppression of CR diffusion can explain the presence of radio ghosts. We show that 
the partial escape of cosmic rays is mostly confined to the wake of the rising bubbles, 
and speculate that this effect could: (1) account for the excitation of the Ha filaments 
trailing behind the bubbles in the Perseus cluster, (2) inject entropy into the metal en- 
riched material being lifted by the bubbles and, thus, help to displace it permanently 
from the cluster center and (3) produce observable 7-rays via the interaction of the 
diffusing cosmic rays with the thermal intracluster medium (ICM). 

Key words: ICM: outflows - MHD - magnetic fields - AGN: clusters of galaxies - 
cosmic rays 



1 INTRODUCTION 

The discovery of X-ray cavities in the ICM that coincide 
with radio emission (e.g., Blanton et al. 2001, Boehringer et 
al. 2003) , sound waves and weak shocks (Fabian et al. 2003a) 
provided clear observational support for the idea that AGN 
are the main heating source that could offset radiative cool- 
ing and prevent "cooling catastrophes" in clusters of galax- 
ies. There is growing consensus that overall heating and cool- 
ing balance may indeed be achieved. This is based on the 
analysis of statistically significant samples of ICM bubbles 
(Birzan et al. 2006) as well as other arguments. However, the 
actual mechanism by means of which the mechanical energy 
of AGN is delivered to the thermal ICM remains elusive. 

The bubbles that supply energy to the ICM are known 
to be filled with magnetized, relativistic, radio-emitting, 
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non-thermal plasma and, thus, the diffusion of cosmic rays 
(CR) from magnetic bubbles and their interaction with the 
ICM is one of the key physical mechanisms that needs to be 
explored in the context of AGN feedback in clusters. These 
processes are of high astrophysical relevance to a number of 
long-standing problems. If the diffusion rates are relatively 
high, cosmic ray diffusion may explain the absence of strong 
"cooling flows", the origin of cluster radio halos and, as we 
suggest in this paper, the origin of the excitation of Ha fila- 
ments in clusters □ and the metal distribution in clusters. If, 
on the other hand, CR remain largely confined to the bub- 
bles as they rise buoyantly, they may explain the presence of 
radio relic sources. Such electrons will be protected against 
Coulomb losses and are thus able to produce a significant 
non-thermal Comptonisation signature of the CMB. 

1 A similar possibility was independenty suggested by Sanders & 
Fabian (2007) 
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The effects of CR diffusion in clusters have been pre- 
viously studied by, e.g., Mathews & Brighenti (2007), who 
used a one-dimansional model to put upper limits on the 
effective diffusivity rate of cosmic rays from AGN-inflated 
bubbles. These effects were also addressed by Jones et al. 
(1999) and Tregillis et al. (2004), who performed 3D mag- 
netohydrodynamical (MHD) simulations focusing on a dif- 
ferent regime (high Mach number jets). Cosmic ray diffu- 
sion was also incorporated recently into cosmological struc- 
ture formation simulations by Miniati et al. (2001) (see also 
Miniati 2001 for the description of the COSMOCR code) 
and by Enfilin et al. (2007) and by Jubelgas et al. (2006) 
using the GADGET code. 

In this paper we consider three-dimensional MHD sim- 
ulations of AGN-inflated bubbles in their buoyant stage. We 
go beyond previous treatments of this feedback mechanism 
and include the anisotropy of cosmic ray diffusion introduced 
by magnetic fields in the bubbles. In the next section we de- 
scribe the methods employed to simulate anisotropic diffu- 
sion. In section 3 we present our results and in the last sec- 
tion we discuss the astrophysical consequences of our find- 
ings. 



2 SIMULATION DETAILS 

2.1 The code 

The simulations were performed with the PENCIL code 
(Dobler et al. 2003, Haugen et al. 2003, Brandenburg et 
al. 2004, Haugen et al. 2004). PENCIL is a highly accu- 
rate grid code that is sixth order in space and third order 
in time. It is particularly suited for compressible turbulent 
MHD flows. Magnetic fields are implemented in terms of a 
vector potential so the field remains solenoidal throughout 
the simulation. The code is memory- efficient, uses Message- 
Passing Interface and is highly parallel. 

2.2 Initial conditions 

The details of the numerical setup pertaining to the ther- 
modynamical state of the cluster atmosphere and mag- 
netic fields were described in detail in Ruszkowski et al. 
(2007) (hereafter R07). Here we summarize only the es- 
sential points and differences. The distribution of magnetic 
fields in all runs presented here corresponds to the "drap- 
ing case" whereby the characteristic lengthscale of magnetic 
field fluctuations exceeds the bubble size. The magnetic pres- 
sure is much smaller than the gas pressure (typical plasma 
/3 parameter is of order 20; see R07 for the discussion of 
this choice). Even though such a magnetic field is formally 
"dynamically unimportant" , it nevertheless strongly affects 
the dynamical state of the bubbles by preventing their dis- 
ruption via the Rayleigh- Taylor and Kelvin-Helmholz insta- 
bilities, thus making their morphologies consistent with X- 
ray observations. Bubbles were initially in total (i.e., CR 
and thermal gas) pressure equilibrium with the ambient 
ICM. For numerical reasons we chose the CR pressure to 
be such that its volume-averaged value was comparable to 
the thermal pressure. More specifically, the cosmic ray pres- 
sure fraction /cr was distributed inside the bubble accord- 
ing to fcn(r) = (1 - / g ascos(0.57rr/ri,), where / gas = 0.1 is 



the pressure fraction due to the thermal gas at the centre 
of the bubble, rb = 12.5 kpc is the bubble radius and r is 
the distance from the bubble centre. We assumed that the 
cosmic ray contribution to pressure outside the bubble is 
negligible. The bubbles were made buoyantly unstable. At 
each point within the bubble volume, the gas density was 
pgas = Picm(1— fcn)/x and the gas temperature T = Ticm^;, 
where x — 10. The temperature of the ambient ICM was 
constant and equal to lOkeV (note that high values of tem- 
perature lower the effective viscosity of the code as compared 
to the Braginskii value). A constant value for the ambient 
temperature also makes it easlier to identify the volume oc- 
cupied by the bubble (i.e., there is no need to introduce a 
dynamically passive fluid to define the bubble boundaries; 
we use this fact in the analysis of our simulations below). 

The resolution of our simulations was 200 3 . We used 
the following code units in this work: code time unit was 
3.3 x 10 6 yr, code length unit was 1 kpc and density 
unit 10~ 24 g cm' 3 . In these units, the diffusivity K = 
32.4(iCo/3 x 10 30 cm 2 s _1 ) code units, where K denotes the 
diffusivity in cgs units. 

2.3 Cosmic ray diffusion 

Cosmic rays are incorporated into the code via a two-fluid 
model. The gas and cosmic ray plasmas have adiabatic in- 
dices equal to 7 = 5/3 and 7 cr = 4/3, respectively. Cosmic 
rays are coupled to the gas and act on it via their pressure 
gradient. 

Diffusion of relativistic particles in PENCIL takes into 
account the effects of anisotropy introduced by magnetic 
fields (diffusion in the direction perpendicular to the field 
lines is slower than along the field lines). 

The standard diffusion equation violates causality by 
allowing for infinite particle propagation speeds. The diffu- 
sion method implemented in the PENCIL code (Snodin et al. 
2006) improves on the standard Fick diffusion law and takes 
into account the requirement that particles should propagate 
at finite speeds. This is achieved by expanding the Boltz- 
mann equation in the relevant smallness parameter and re- 
taining terms up to second order (see Gambosi et al. 1993 for 
the discussion of charged particle transport). This treatment 
is an extension of the method of Hanasz & Lesch (2003), who 
implemented CR diffusion in the ZEUS code, that included 
anisotropic diffusion but did not consider departures from 
the standard diffusion law. The relevant smallness parameter 
in the expansion of the Boltzman equation is the Strouhal 
number St = (K yr) 1 / 2 //, where (K\\ is the diffusion coef- 
ficient along the magnetic field lines, I is the characteristic 
lengthscale of the initial magnetic structure and r is the 
damping time of the non-Fickian contribution to diffusion; 
Landau & Lifshitz 1987) . This approach leads to the general- 
ization of the diffusion equation to the "telegraph" equation 
that reduces to the standard diffusion equation in the limit 
of St — > but becomes a wave equation for St — > 00. A non- 
Fickian diffusion model has also been derived in the context 
of turbulent diffusion of passive scalars by Blackman & Field 
(2003) and confirmed in numerical experiments (Branden- 
burg, Kapyla & Mohammed 2004). 

A very rough estimate of the diffusion coefficient can 
be obtained from a dimensional argument based on the 
radius of the inner X-ray cavities in the Perseus cluster 
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r ~ 5r5 kpc and the time it takes to inflate the bubble 
t ~ 10 7 i 7 yr Mathews & Brighenti 2007). These values yield 
^ < 7.5 x 10 29 r|/t7 but the estimate is somewhat uncer- 
tain. This value is consistent with the theoretical estimates 
of Enfilin (2003) /ty ~ 2 x lO^S^r^-Br^^cmV 1 , 
where Uio is the energy of radio emitting electrons in units 
of 10 GeV, Bi is the magnetic field in /xGauss, and 8 is a 
parameter of order unity or less that depends on the power 
spectrum of magnetic field fluctuations and the ratio of par- 
ticle scattering frequency to gyro-frequency. Assuming that 
r ~ r/vA, where va is the Alfven speed, r is of order unity 
in code units for the typical parameters in our simulations. 
As r is essentially a free parameter of this model we adopt 
this value here. Guided by the above estimates for the dif- 
fusion coefficient and the damping time, we estimate the 
Strouhal number to be St ~ 0.5(tK ]} / 'K) 1/2 (r /n)' 1 , where 
rb — 12.5 is the bubble radius in code units and K — 32.4 
is the reference diffusivity in code units (in cgs units it is 
3 x 10 cm s~ ). The code requires some magnetic diffusiv- 
ity, rj, viscosity, v, and an isotropic contribution to cosmic 
ray diffusivity, Ki BO , to ensure numerical numerical stability. 
These coefficients can be lowered at the expense of increas- 
ing numerical resolution. Except where explicitly stated, in 
our simulations we adopted rj = v — K 1BO = 0.07. As we 
intend to study the effect of the anisotropy of cosmic ray 
diffusion, we want to keep the ratio of K\\/Ki BO as large as 
possible. We also want to be conservative in our assessment 
of the importance of the diffusion and, therefore, consider 
high values of diffusivity. We stress that for high diffusivi- 
ties, simulations would not be feasible if not for the fact that 
the code uses a non-Fickian approach, i.e., effectively limits 
the CR propagation speed and prevents the timestep from 
decreasing to unacceptably low values. 



3 RESULTS 

Figure 1 shows snapshots from our simulations in the high 
diffusivity regime. The first column shows the natural loga- 
rithm of thermal gas density in the cluster, while the second 
and third columns show the logarithm of the CR energy 
density for two different assumptions about the isotropy of 
diffusion. The second column corresponds to the case where 
diffusion in the direction perpendicular to the local magnetic 
field is suppressed (i.e., K± = 0, K\\ = K), while the last 
column is for the case where the diffusion coefficients of cos- 
mic rays in the direction parallel and perpendicular to the 
local direction of magnetic field lines are identical. As the 
density distribution looks quite similar in both the isotopic 
and anisotropic cases, we show only the density map for the 
anisotropic diffusion case. The reason for this similarity is 
that the pressure support in cosmic rays is comparable to 
that in the thermal component (for numerical stability rea- 
sons) and, thus, even when all cosmic rays diffuse out of the 
bubble, its radius diminishes only by a modest amount of 
~ 1 - OM-W^ ~ 10%. 

The first column in Figure 1 demonstrates that the 
"draping effect" of the magnetic field prevents the disrup- 
tion of the bubble as it rises buoyantly in the cluster atmo- 
sphere. We also note that the magnetic field does not show 
any signs of decay on the timescales considered here, either 
due to numerical effects or to natural relaxation of magnetic 



fields. While the bubble maintains its integrity, the cosmic 
rays gradually diffuse out of it. However, this process is not 
isotropic. In the second column we see that the spatial distri- 
bution of cosmic rays initially tends to follow the distortions 
in the magnetic field generated by the upward motion of the 
bubble (lower panel, second column). Magnetic fields and 
cosmic rays follow the velocity vortex pattern of the rising 
bubble. At later times (upper row), the magnetic fields in- 
side the bubble "open up". This, however, takes place only 
in the wake of the rising bubble. The working surface of 
the bubble is efficiently protected from disruption by a thin 
layer of ordered magnetic fields (see R07 for more detailed 
discussion of how the field topology changes due to bubble 
motion). These ordered fields not only prevent instability 
from developing but also stop cosmic rays from escaping the 
bubble. On the other hand, magnetic fields stretch out be- 
hind the bubble and become approximately parallel to the 
direction of the motion. Cosmic rays can therefore diffuse 
more easily in the direction of gravity. This effect is clearly 
seen in the upper row of the second column. In the last 
column of Figure 1 we show what happens when diffusion is 
isotropic (diffusion rates in the directions perpendicular and 
parallel to the magnetic field lines are identical and equal to 
the parallel diffusion rate shown in the second column). In 
this case, the diffusion time is much shorter than the buoy- 
ancy time and the pattern of energy density in cosmic rays 
is markedly different from that corresponding to anisotropic 
diffusion. In the initial stages of the evolution of CR energy 
density a wave- like pattern is clearly seen - the diffusion sig- 
nal propagates at finite speed. 

Figure 2 shows the same quantities as Figure 1 but for 
the case where the fiducial magnitude of the diffusion coeffi- 
cient K is reduced by a factor of ten. The viscosity coefficient 
v was increased to 0.2 to prevent numerical instabilities for 
the runs presented in this figure. The arrangement of panels 
in this figure is the same as in Figure 1. The density evo- 
lution is practically identical. The comparison between the 
second column with its equivalent in Figure 1 shows that, as 
expected, cosmic rays are now better confined to the bub- 
bles (compare the pattern shapes and the CR intensities). 
The third column shows that although diffusion is formally 
isotropic (K± = K\\), the diffusion rate is slow enough for 
the fluid to advect cosmic rays further away from the cluster 
centre. 

In Figure 3 we quantify the above results in terms of the 
integrated cosmic ray energy normalized to the total initial 
energy inside the bubbles as a function of time. Dashed lines 
correspond to the high diffusivity case (cf. Figure 1) and 
solid lines to low CR diffusivity (cf. Figure 2). In each of 
these two cases, the lower curve corresponds to the isotropic 
case and the upper one to the case of anisotropic diffusion. 
These curves clearly demonstrate that the anisotropy effects 
introduced by the topology of the magnetic field are cru- 
cial in determining the escape fractions of cosmic rays from 
AGN bubbles. Anisotropic diffusion can lead to a substan- 
tial fraction of cosmic rays being retained inside the bubbles 
on timescales of order the buoyancy time. Any diffusion that 
takes place in such cases is confined to the wake of the bub- 
ble. We now move on to discuss the implications of these 
results for the physics of the intracluster medium. 
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4 ASTROPHYSICAL CONSEQUENCES 

The above results have the following astrophysical conse- 
quences: 

• Radio relics 

Our simulations show that while diffusion out of the bubbles 
does take place, a significant fraction of the original energy 
in cosmic rays can be confined to the buoyant bubbles, 
provided that the perpendicular diffusion coefficient is 
not too large. Lower values of diffusivity than those 
considered here would only strengthen this conclusion. 
If this is the case, then cosmic rays can be efficiently 
screened from Coulomb losses due to their interaction with 
the ambient medium. Such energetic plasma confined in 
the bubbles, when subject to passing waves, (either of 
AGN origin or due to infalling substructure clumps) could 
be re-energized and lead to detectable radio emission, 
thus explaining cluster "radio relic" sources (Enfilin & 
Briiggen 2002). The fraction of cosmic rays that does 
escape radio cocoons will pollute the ICM and contribute to 
the so called "radio haloes" Enfilin & Gopal-Krishna (2001). 

• Abundance profiles in clusters of galaxies 

X-ray observations (Schmidt et al. 2002, De Grandi et al. 
2004, Bohringer et al. 2004; in Virgo, Perseus, Centaurus 
and Abell 1795 clusters reveal that clusters without cool 
cores show a spatially uniform distribution of metals while 
cooling flow clusters reveal strong centrally peaked metal- 
licity profiles. However, the observed metallicity profiles 
appear much broader than the stellar light distribution 
of the brightest cluster galaxies (BCGs). One possible 
solution to this puzzle is that AGN help to uplift the gas 
enriched with metals (e.g., Briiggen 2002, Roediger et al. 
2007). In the absence of any microscopic energy transfer 
between the interior of the bubble filled with relativistic 
plasma and the displaced metal-enriched thermal gas, the 
latter will tend to sink back toward the cluster center. This 
problem is closely related to the problem of the develop- 
ment of hydrodynamical instabilities on the bubble-ICM 
interface. In the absence of any mechanism to suppress 
such instabilities (such as viscosity or magnetic fields), the 
rising bubble will be quickly shredded and the non-thermal 
bubble material will mix with the ICM, adding entropy 
to the metal-enriched gas and preventing it from sinking 
toward the cluster center. This is what happens in pure 
hydrodynamical simulations. In a realistic situation, mag- 
netic fields may help to prevent such instabilities and the 
entropy gain of the metal-enriched gas must have a different 
origin. Here we suggest that the diffusion of cosmic rays out 
of the bubbles and their subsequent thermalization may 
provide the necessary energy injection. Depending on the 
energy of radiating electrons and the density of the ICM, 
either Coulomb interaction of CR electrons with the ICM 
(^cooT' ~2x 10 s 727^^2 y r ) or synchrotron/inverse Compton 
energy losses will dominate {t' C ao\ ~ 2 x 10 8 7j 1 [l + ub/uic] 
yr), where n_2 is the electron density in units of 10 -2 
cm~ 3 , 72,4 is the relativistic gamma factor in units of 10 2 
and 10 4 , respectively, and ub,ic are the energy densities 
in magnetic fields and CMB photons, respectively. The 
thermalization of CR electrons by Coulomb interactions 



may then preferentially increase the entropy of the metal 
enriched material being pulled up in the cluster atmosphere 
by buoyantly rising bubbles. This is because cosmic rays 
tend to preferentially diffuse in the direction opposite to 
their motion, i.e., their thermalization will mostly take place 
in the metal enriched ICM dragged behind the bubbles (at 
least for the magnetic field parameters considered here). 
The observations of high abundance ridges displaced by 
fossil bubbles reported by Sanders, Fabian & Dunn (2005) 
are generally consistent with this idea. 

• Excitation of Ha line in the filaments trailing behind 
bubbles 

Emission-line nebulae are commonly found surround- 
ing massive galaxies in the centres of cool cluster cores 
(Crawford et al. 1999). The best examples of extended 
Ha-emitting filaments pointing from the cluster centre to 
the buoyant bubbles are those in the Perseus cluster (Hatch 
et al. 2006). Several mechanism have been proposed (and 
ruled out) to explain the excitation mechanism of this 
line. One possibility is that these nebulae are excited by 
stellar UV, but there is no spatial correlation between the 
filaments and stellar clusters (Hatch et al. 2006). Another 
mechanism could be excitation by X-rays from the ICM 
but the filaments are up to a hundred times less luminous 
in X-rays than in the UV (Fabian et al. 2003b). Finally, it 
has been proposed that heat conduction from the ICM to 
the colder filaments could provide line excitation (Donahue 
et al. 2000). However, thermal conductivity would also lead 
to the evaporation of these filaments (Nipoti & Binney 
2004). Here we suggest the possibility that cosmic rays, 
preferentially diffusing along the field lines trailing behind 
the rising bubbles, could provide the excitation mechanism 
for the filaments that are also located in the bubble wakes. 
We note that the same magnetic fields that act like "wires" 
conducting cosmic rays could also prevent the filaments 
from evaporating via thermal conductivity by locally 
suppressing it. As a side note we add that while the ICM 
may be viscous, the laminar appearance of the filaments 
does not necessarily imply that the flow is viscous. The 
ordered and amplified magnetic fields trailing behind the 
bubble may prevent the destruction of the filaments by 
turbulent motions. 

• j -ray signatures 

It has recently been proposed that extended radio features 
of radio galaxies (Sambruna et al. 2007; in case of Fornax 
A) and the interaction of cosmic rays from AGN cavities 
with the ICM (Hinton & Domainko 2007; Hydra A) could 
be detectable with GLAST. Moreover, it has been pointed 
out that the interaction of such CR with the ICM could lead 
to significant annihilation line emission (Furlanetto & Loeb 
2002, Griffiths 2005). We point out that our simulations add 
more weight to this idea and that a preferential diffusion of 
CR in the wakes of rising bubbles could potentially lead to 
observable 7-ray signatures. 

We note also that these results emphasize the role of 
non-ideal hydrodynamical effects in the studies of AGN feed- 
back in general. While we argued here that cosmic ray diffu- 
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Figure 3. Total energy in cosmic rays inside buoyant bubbles 
normalized to the total initial energy of cosmic rays as a function 
of time. Dashed curves are for the high diffusivity case (cf. Figure 
1) and the solid ones are for low diffusivity (cf. Figure 2). In each 
of these two cases, lower curves are for isotropic diffusion and the 
upper ones for suppressed diffusion in the direction perpendicu- 
lar to the local orientation of magnetic fields. See text for more 
details. 

sivity may be one of the key effects that needs to be included, 
other transport process such as conduction (Ruszkowski & 
Begelman 2002, Bruggen et al. 2005, Fabian et al. 2005) or 
viscosity (Fabian et al. 2003a, Ruszkowski et al. 2004a,b, 
Reynolds et al. 2005, Sijacki & Springel 2006) may also play 
an important role. Finally, if the "pollution" of the intra- 
cluster medium by cosmic rays diffusing out of AGN bubbles 
is significant, then other mechanisms such as the generaliza- 
tion of the magnet othermal instability (Balbus 2004, Parrish 
& Stone 2005) by Chandran (2005) may also play a role. We 
note that recent X-ray observations with Chandra and XMM 
Newton show evidence for nonthermal particle populations 
in the ICM (Sanders et al. 2006, Werner et al. 2007). 
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natural logarithm of density, while the second and third columns present the logarithm of the energy density of cosmic rays for vanishing 
diffusion in the direction perpendicular to the local magnetic field (Kj_ = 0, Ku = K; middle column) and for the case where the diffusion 
of cosmic rays in the direction parallel and perpendicular to the local direction of magnetic field lines are identical (K± = Ku = K). 
Lower row corresponds to earlier times t = 3, upper to t = 25 (one code time unit is 3.3 X 10 6 years). In the case of density, the color 
bar extends from -7.0 to -2.0 in code units. All the remaining columns are for the range [-3.0, 0.0]. See text for more details. 
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Figure 2. Same as Figure 1 but for the fiducial value of diffusivity K decreased by a factor of 10. 
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